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ABSTRACT 

Gravitational amplification of Poisson noise in stellar systems is important on large 
scales. For example, it increases the dipole noise power by roughly a factor of six and 
the quadrupole noise by 50% for a King model profile. The dipole noise is amplified by a 
factor of fifteen for the core-free Hernquist model. The predictions are computed using 
the dressed-particle formalism of Rostoker & Rosenbluth (1960) and are demonstrated 
by n-body simulation. 

This result implies that a collisionless n-body simulation is impossible; The fluc- 
tuation noise which causes relaxation is an intrinic part of self gravity. In other words, 
eliminating two-body relaxation does not eliminate relaxation altogether. 

Applied to dark matter halos of disk galaxies, particle numbers of at least I0 6 will 
be necessary to suppress this noise at a level that does not dominate or significantly 
affect the disk response. Conversely, halos are most likely far from phase-mixed equi- 
librium and the resulting noise spectrum may seed or excite observed structure such 
as warps, spiral arms and bars. For example, discreteness noise in the halo, similar to 
that due to a population of 10 6 M Q black holes can produce observable warping and 
possibly excite or seed other disk structure. 

Key words: gravitation — methods: numerical — methods: statistical — stellar 
dynamics — galaxies: kinematics and dynamics — galaxies: evolution 



1 INTRODUCTION 



Stellar systems are finite number systems by nature. Because 
the characteristic relaxation time is so long in many cases of 
astronomical interest, the continuum limit is a good approxi- 
mation and the collisionless Boltzmann equation (CBE) gov- 
erns the evolution. Because the CBE is difficult to solve 
in any generality, many researchers have turned to n-body 
simulation as a primary source of insight. For a modest 
number of particles, however, the naive n-body problem is 
not a simulation of the CBE but of the collistonal Boltz- 
mann equation. Enormous effort has been applied to ma- 
nipulating the interparticle force to reduce the intrinsic re- 
laxation to produce a near- collisionless solution. The most 
commonly used technique smoothes the particle over a finite 
size region leading to the softened the point mass potential: 
Gfx/r — > Gfi/V r 2 + a 2 . As long as the smoothing size is not 
larger than the mean interparticle spacing, intuitively, there 
should be no significant change in the results. 

Unfortunately, two-body interactions are only part of 
the story. Poisson fluctuations excite structure at all scales 
in the system. Many simulations are optimized to resolve 
large-scale features but the relaxation is enhanced by large- 
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scale collective excitations on these same scales (Weinberg 
1993). One can think of this excitation as the projection 
Poisson noise on the modal spectrum of the stellar system. 
This is the same modal spectrum responsible for producing a 
response by a perturber of interest, such as a galaxy reacting 
to orbiting or passing companion. Therefore, a collisionless 
solution is physically impossible without irrevocably chang- 
ing the dynamics of the system under study. In other words, 
one suppresses the global part of the relaxation at the risk 
of throwing out the physics responsible for the evolution one 
is studying. There is no other recourse but large numbers of 
particles. 

For example, this work was motivated by n-body ex- 
periments with a thin disk in a live halo. The observed fluc- 
tuations vertical in force at the disk plane were larger than 
predicted for Poisson fluctuations for n-body simulations of 
halos with 10 particles using the SCF scheme (e.g. Clutton- 
Brock 1972, 1973, Hernquist & Ostriker 1992). However, 
real galactic halos are certainly not smooth and contain gas 
clouds, star clusters, dwarf galaxies, stellar streams, and pos- 
sibly as of yet undetected massive objects such as 1O 6 M0 
black holes (e.g. Lacey & Ostriker 1985). All of these can 
contribute to correlated fluctuations at large scales leading 
to warped disks and other possibly observable distortions 
(see Weinberg 1997). 

In this paper, I describe the expected amplitude of 
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noise-generated fluctuations in a spherical equilibrium stel- 
lar system including self-gravity. The main result is the 
power spectrum of fluctuations generated by long-range cor- 
relations of particles moving their own gravitational field. 
This is computed using the polarization cloud method devel- 
oped by Rostoker & Rosenbluth (1960) for plasma physics. 
The same approach can be used for any regular system (see 
Nelson & Tremaine 1997 for a general discussion). The power 
at very small scales will be Poisson, but at large scales, it 
will be modified by the global gravitational response. The 
basic results developed in §^ show that this might have been 
predicted a priori and is applied to a few standard spherical 
models in §ra and the predictions corroborated by n-body 
simulation. We conclude in 



2 DERIVATION 
2.1 Basic approach 

Since we are concerned with perturbations to an equilibrium, 
we may solve the linearized collisionless Boltzmann equation 
(LCBE): 



dfi dfi 0H o df gflj _ Q 
dt dw dl dl dw 



(1) 



where the subscripts o and 1 denote the background and 
first-order perturbed quantities and H is the Hamiltonian. 
The LCBE has been written in action-angle variables. 

I will solve the LCBE by a Laplace transform in time 
and a Fourier transform in angle variables (Tremaine & 
Weinberg 1984, Weinberg 1989). Let us denote the Laplace 
transform of some quantity q(t) by q(s) and Fourier trans- 
form of q by qi (I). The vector quantity 1 is the triple of 
'quantum' numbers defining the discrete Fourier series in 
angle variables. Recalling that the frequencies corresponding 
to the angle variables are CI = dH /dl, the Fourier-Laplace 
transform of the LCBE (eq. |l|) becomes 



■A + a-nA-a~#u = o. 

ol 

Solving for f\ yields 



A = ii . 



df Hu 
dl a + il-n' 



(2) 



(3) 



2.2 Dressed particles 

We will use the solution in equation to derive the re- 
sponse of the continuum stellar system to point particles. 
The effect of a point particle on the system is then the com- 
bined effect of the potential due to point particle and its 
response. This has been called a dressed point particle by 
Rostoker & Rosenbluth (1960) who first derived the proper- 
ties of a plasma of dressed particles. 

Following Weinberg (1989), I will expand the pertur- 
bation in a biorthogonal series whose basis is constructed 
from eigenfunctions of the Laplacian. The potential or den- 
sity, then, trivially follow from the expansion coefficients and 
the potential-density pairs, Ui,di, solve the Poisson equa- 
tion explicitly. I will deviate from traditional notation and 
define the true space density corresponding to m to be 



di — —di/AnG. The biorthogonal expansion, then, is in po- 
tential and — 47rG times the density. This leads to the con- 
venient orthogonality condition 

drr 2 Ui U *d l j n = Sij. (4) 

The upper limit of the integral in equation (0) may be infi- 
nite, depending on the pairs. I will set G = 1 throughout. 

Using this, the spherical harmonic expansion coeffi- 
cients for a point mass of mass fa — M/N at r = r (t) 
is 



= -A-Kfh \ d 3 r'Y l l 



, :nK 9' ,4>')u- n *{r')8 3 {y' - r 



-A-KmY* m {0 o , <j> )uj n *(r )- 



(5) 



The gravitational potential of the point mass, the perturbing 
Hamiltonian, is then 



(6) 



l,m j 



2.3 Expansion in action-angle series 

To solve the LBCE using equation (^), we now need to ex- 
pand functions of the form Yi m {6, <f>)g(r) in action-angle vari- 
ables. Specifically, for a spherical system, it is easy to write 
down the description of a particle orbit in the orbital plane. 
From this, we may derive the general harmonic expansion 
following the technique presented in Tremaine & Weinberg 
(1984). This yields 



Y lm {e^)g{r 



= 8„ 



oo 



I 



^ jm+i 2e il.w x 

t = — oo lo = — l 



F„ 2 (7r/2,0)rl 2m (/3)W# 2m (I 



(7) 



where (3 is the elevation of the orbital plane defined by the 
actions I, W is 

Wll m {l) = zl dwicoBilim+hfimMrivn)) (8) 



with f(wi) = ip — W2, and r\ ((5) is rotation matrix for 
spherical harmonics (e.g. Edmonds 1960). 

2.4 The Laplace transform of Hu 



Putting the results of §§2.1-2.3 together, we can derive the 
Laplace transform of Hu and using equations (^) — (^) to 
get the distribution function of the dressed particle. 

First, for a particular spherical harmonic, the action- 
angle transform of Hi is 



Hi 



(2*) 



j3 — zl-w X ^ i I'm 

awe \ i, 



/ J b 3 m Y lm (9, 4>)u- n (r) 



= ^y ii2 (V2,0)r| 2m (/3X^(I)fefW 



(9) 



where &j m (t) is a some general well-behaved function of time. 



The quantity Wjj^QO is shorthand for W with g(r) = 
«j- m (r). Because the motion of a star on a regular orbit is 
quasi-periodic, we will consider terms b l J m (t) with pure sinu- 
soidal dependence: fej m (f) = b 1 ™ 1 exp(iujt) . Substituting this 
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into equation (g) and Laplace transforming gives the desired 
result: 



where the Laplace transform of bj(t) = b ™ exp(iujt) is 



(11) 



2.5 The response of the dressed particle 



We can compute the density and potential response cor- 
responding to the dressed particle by integrating the per- 
turbed distribution function from equation (||) over veloci- 
ties: 



p\{s) 



d 3 vf\. 



(12) 



The Laplace transformed expansion coefficients of the 
biorthogonal basis are then 



a = -4tt / d rY t 



i\ m *(r) I d 3 vf_ 



(13) 



This computation is easily performed by noting that the Ja- 
cobian of the canonical transform of d 3 rd 3 v is unity and we 
are free to choose any set. This procedure is the motivation 
behind the biorthogonal expansion. If we choose d 3 rd 3 v = 
d 3 Id 3 w and use equation ([?]), we can do the angle integration 
trivially. Then, noting that d 3 I = dEdJJd(cos (3)/fli(E, J), 
we can do the integral in (3 using the orthogonality of the 
rotation matrices: 



2l + l S ^ Sl3 



(14) 



(Edmonds 1960). We get: 
a\ m = -47r(27r) 



21 + 



dEdJJ ., df 



Ui{E,J) 81 



= — / dse st V- llm {s)n lm (s)-h lm — — . 

ZTli / s — IU1 

</ c — zoo 

(18) 

Assuming that the background is stable with no oscillatory 
modes, T>~ llm (s) has no poles in the half plane 5ft (s) > 0. 
The final term gives a pole at s = iu). Although there may be 
poles in the half plane 5R(s) < 0, these will vanish for t — > oo 
relative to the pure imaginary contribution. To perform the 
integral, one may take the s integration into the phase-space 
integral for the elements of TZ lm (s). In addition to V~ llm (s), 
the s-dependence is in two simple poles and one finds an 
integral of the form 



1 

2ri, 



7 st -r>— 1 lm / \ 

dse V (s) 



TX 



C — ZOO 

1 lm f • A Aujt 

(iui)e 



V 



s + il ■ £1(1) s — iu) 
1!m (-il-n(I))e- il ' n(1)( 



■il - n(i) 



(19) 



For large values of t, this expression oscillates rapidly and we 
may extract the dominant coherent contribution. There are 
two cases: without and with a resonance in the phase space. 
The existence of a resonance in phase space is defined by 
u + 1 • £1(1) = for /(I) 7^ 0. For the non-resonant case, 
the integrand has no singularity and we can consider each 
term separately. The first term yields a contribution in phase 
with the perturbation while the second term in the brackets 
oscillates incoherently and makes no net contribution. The 
second term, therefore, can be ignored. For the resonant 
case, the contribution at large t has a sharp peak about ui + 
1 ■ £1(1) = as t — > oo. Expanding T)~ llm (iu) about iuj and 
retaining only dominant terms, one finds the contribution 
near the resonance is 



r i T1 -lim/. n i(ui — l-0(I))t/2 „ ' 

] » T> (icu)e K K " 1 2 ■ 



j+i-n(i) 



or 

lim [ ] = V- llm (iu)e iwt 2TrS(v + l-£l(l)). (20) 

t — >oo 



s + il-£l 



\Y ll2 (tt/2, 0)|X,V* (W a ' TO W 



(15) 

Equation @ has the form a' m = -R lm h lm and de- 
scribes the response of the stellar system, a lm , to the per- 
turbation, b lm . The self-gravitating response, then, is the 
solution to a lm (s) = TZ lm (s) ■ [a lm (s) + h lm (s)] : 



a im (s) = V- llm {s)1l lm {s)-h lm {s) 
where 

v tm (s) = i-n lm {s). 



(16) 
(17) 



If b = 0, the equation for the response becomes an eigen- 
value problem. The eigenvalues s are the zeros of the disper- 
sion relation detV lrn (s) — 0. For the general stable spherical 
stellar system, there will be no modes for 5R(s) > and only 
in special cases with restricted phase space does one find 
oscillatory modes, 5ft(s) = 

To evaluate the coefficients as a function of time, we 
perform the inverse Laplace transform, deforming the con- 
tour to the 5R(s) — > — oo: 



a im (t) = _!_ / dse st V- llm (s)n lm (s)-h l 

J c — zoo 



We will adopt the latter asymptotic form here and see in 
the final computation that the exp(iurf) will cancel leaving 
only the delta functions. Putting both cases together yields 
a simple expression 



a-i (t) = 



2_^M iv (u)b v , 



(21) 



-47r(27r) 3 



21 

dEdJJ 



1 ^ 



dfo 



[iu) 



i 



£li(E, J) 81 

f„ 2 (V2,o)| 2 </ i ;(iX;ji)x 
i 



V 



(w + i-n(i)) 



+ 2ni5 (w + l- n(I)) 



(22) 



To simplify notation, we have explicitly noted in equa- 
tions @ and ( |22| ) that the solution takes the the form 
a lm (t) = exp(iujt)M Lm {Lu) ■ h lm . The integrals in the ma- 
trix elements of T> lm may be analytically continued using 
the Landau prescription (e.g. Krall & Trivelpiece 1973) af- 
ter a conformal mapping of the discrete interval in E. Notice 
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that this expression takes both resonant and non-resonant 
cases into account; without a resonance, the delta function 
does not contribute and principal value of a non-singular 
integrand is the integrand itself. One recovers the non-self- 
gravitating but global response by setting T>^ — <5y. 



2.6 Energy in fluctuations 

We will now use equation (^l|) to evaluate the fluctuation 
energy at different spatial scales assuming that individual 
particles are uncorrelated. The particle wakes do in fact give 
rise to correlations but this is of higher order in 1/N in the 
BBGKY expansion (cf. Gilbert 1969) than the lowest-order 
effect we will consider here. 

This leaves us with individual particles reacting coher- 
ently to the effect of their own wakes. Because the particles 
are uncorrelated, the number density of particles at I ,w 
at time and at I', w' at time t is 

P(I o ,w o ,0;l',w',t) = P(I ,w )5(l'-I ) X 
5(w - w' + U[I )t) 



Gathering terms, this can be simplified as follows: 



(23) 

where V(l , w ) is the equilibrium particle distribution with 
N= I d 3 I d 3 w V(I ,w ). (24) 



Direct substitution demonstrates that equation (p3|) solves 
the Liouville equation with the initial condition I = I D and 
w' = w at t = 0. Similarly, integrating equation (^3|) over 
all coordinates gives N. 

For a given harmonic Im, the fluctuation energy is then 

~<<&ipr) = lJd 3 r(^2aT f Y i : n (e, ( p)u l r'(r)x 

i 

7~ E a * my i™( 61 ' <t>)d l ™{r)) 



1 / 1 \ / lm* lm\ 



(25) 



where the expectation value of some quantity H is defined 
by 



(S) = / d 3 I o d 3 w o d 3 l'd 3 w'T{Io,Wo,0-,l',w',t)E 

= N J d 3 I d 3 w d 3 l'd 3 w'f (l )S{l'-l )x 
5(w — w' + fi(I )t) H. 

(26) 

Applying equations (^l|) and ( pj| ) to equation ( [iB] ) gives 

2/+ i)EE x 



1 / ^ i 4-Km f . q f 2 
2<*iPi> = ~ { ] 



1 i{iv 



dEdJJ 



Y lh (7T/2,0)\ 2 f o (E,J) x 



ni(E, j)' 

wit£{i)wii: m {i). 



4irm 



dEdJJ 
^(E, J) 



Y lh (n/2,0)\ 2 f o (E,J) 



(28) 



Note that each term in the fluctuation energy, equation (j28|), 
is negative definite as expected. The contribution for each 
triple in the angle expansion, 1, and each term in the basis 
expansion k may be tabulated separately. 

We may compute the fluctuation energy in the 
absence of gravity by returning to equation ( [iB] ) 
and evaluating (a' m *aj m ) without any dynamics. For 
TV particles, the sample value for a 1 ™ 1 is a™ = 



defined by equation 



u ; j m *(rfc). Using the expectation 
one finds: 



1 , . . 4irm 
2<^pi> = - — 



Y] J d i rp{r)u l i 



^(r)u l r{r). (29) 



(27) 



This is identical to equation ( g8J ) with out the particle dress- 
ing: M l ff = S iu . 



3 EXAMPLES OF FLUCTUATION 
SPECTRUM STANDARD MODELS 

We will look at both King models and core-free Hernquist 
models. First, we apply equation to a King model 

(Wo = 5) for harmonics I = m = 1, 2, 3, 4 and compare with 
a SCF simulation for 100,000 particles and l ma x = 4. Figures 
|l]and| show the fluctuation energy per expansion function 
plotted against the index of expansion functions (cf. Fig. ^|) 
and the cumulative fluctuation energy less than given index. 
One energy unit is equal to the total gravitational potential 
energy of the unperturbed sphere and m = 1. There are 
slight systematic differences between the n-body results and 
the predictions but all-in-all, the n-body simulation follows 
the analytic predictions fairly well and confirms the excess 
power at low order due to amplification by self gravity. This 
excess is illustrated in Figure H which shows the total energy 
for harmonic orders / = 1-4 derived from the simulation. 
For 1 = 1, the enhancement is roughly a factor of 6. The 
large magnitude for I = 1 is due to the stochastic excitation 
of a weakly-damped mode. For I = 2, the enhancement is 
roughly a factor of 1.5. For I > 4 the power enhancement 
due to self gravity is negligible. The size of coherent struc- 
tures decrease with increasing harmonic order and therefore 
higher harmonics have power at smaller and smaller scales 
for which self gravity is less important. The index at peak 
power increases with I as expected. 

The case for the core-free Hernquist model is shown 
in Figure ^| for harmonics I = m = 1,2. The profiles are 
more sharply peaked about n = 1 because the expansion 
functions well matched to the model profile (Hernquist & 
Ostriker 1992). Especially for the I = m = 1 harmonic, 
the agreement between the expansion and the simulation is 
better in this case. This is probably due to the choice of 
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Figure 1. Top row: Fluctuation energy with orthogonal func- 
tion index in units of background potential energy. Bottom row: 
Cumulative fluctuation energy with orthogonal function indices 
less than n. The left and right columns show the I = m = 1 
and I = m = 2 harmonics respectively. The spatial profile of 
the potential functions is shown in Fig. [ij for I = 1,2. The solid 
and dashed line follows from the analytic calculation for the self- 
gravitating response and Poisson fluctuations alone, respectively. 
The open circles are computed from the time-averaged expansion 
coefficients in the n-body simulation. 




Figure 2. As in Fig. jljbut for the I = m = 3 (left) and I = m = 4 
(right) harmonics. 



Figure 3. Cumulative fluctuation energy for all harmonics of 
given I (labeled) with orthogonal function indices less than n. 
Shown in units of background potential energy as in Fig. |l]. 

expansion functions. For the I — m — 2, the power appears 
systematically high at larger radial order. As for the King 
model, the power at I = 1 has the largest enhancement by 
self-gravity, now by roughly a factor of 15 and this amplitude 
is verified by the simulation. The enhancement at I = 2 is 
similar, roughly a factor of 1.5. 

The root energy indicates the expected magnitude of 
the density or potential fluctuation and can be multiplied 
by yJl/N to estimate the magnitude for an N particle sim- 
ulation. For galaxian disk embedded in a massive halo, an 
large-scale 0.1% distortion in the halo can have interesting 
consequences for disk evolution (cf. Weinberg 1997). In or- 
der to realistically test dynamical hypotheses for disk-halo 
interactions, we need to suppress noise below this level. This 
requires live halos with N > 10 7 particles for S/N> 10. 

In particular, the fluctuating dipole (I = 1) force field 
differentially accelerates and bends the disk, causing thick- 
ening. The quadrupole (I = 2) in a arbitrary orientation 
warps the disk, causing the sort of warp discussed by Wein- 
berg (1995, 1997) but now due to noise rather than a satellite 
wake. Even without self gravity, discreteness noise is suffi- 
cient to require N > 10 6 . 



4 CONCLUSIONS AND DISCUSSION 

This paper considers discreteness noise and its amplification 
in n-body simulations with particular attention a halo's ef- 
fect on an embedded disk. The overall conclusions are as 
follows: 

(i) The SCF n-body algorithm does not suppress Poisson 
fluctuations due to finite number of particles and Poisson 
fluctuations alone can result in significant distortions. 

(ii) The low-order response is amplified by the coherent 
self-gravitating response of the entire system. 



© 1997 RAS, MNRAS 000, §-| 



4 6 8 2 4 

r r 

Figure 4. Potential functions used in fluctuation energy calculation in Fig. [IJ 
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Figure 5. As in Figure hi but for the Hernquist model. 



(iii) As discussed in Weinberg (1994), spherical systems 
have weakly damped modes at harmonics I = m = 1, 2. 
These are excited by noise and enhance the fluctuation 
power. This is clearly seen in Figures ^ and ^[ Detailed agree- 
ment between the linearized solution and simulation reaf- 
firms the importance of these modes to the self-gravitating 
response. 

(iv) For N particles, the I = m = 1(2) harmonics have 
roughly 4/N(1.5/N) times the energy of the background for 
the W = 5 King model and 15/N(1.5/N) for the Hernquist 



model. The amplification does not rely on the existence of 
a core. 

(v) The noise amplitude in an n-body simulation with 
N = 10 s is comparable to the the amplitude required to 
produce a warp (Weinberg 1997). This analysis suggests sim- 
ulations with N ;> 10 7 will be necessary to achieve healthy 
a signal to noise ratio. Conversely, a halo with N = 10 5 
corresponds to masses of roughly 2 to 6 x 10 6 Mq . This is 
comparable to the black hole mass required to produce the 
disk scale height by the Lacey & Ostriker (1985) mechanism. 
Fluctuations at this level are likely to excite bending modes 
in addition to increasing the disk scale height. 

(vi) Although the fluctuations described in Figure [I] in- 
clude self-gravity, self-gravity is a small (< 10%) perturba- 
tion for harmonics I > 4. 

The gravitationally amplified noise described here is un- 
likely to be a significant contribution to the overall evolution 
of a relaxing system, such as a globular cluster. However, the 
noise-excited I = 1 structure may be sufficient to offset the 
cluster halo from its core. This could decrease the lifetime 
of any loss-cone population and subsequently decrease the 
rate of core cooling. 

This calculation investigates the magnitude and nature 
of halo fluctuations assuming that the halo is an equilib- 
rium phased-mixed distribution. However, the outer parts 
of galaxies can barely count 10 dynamical times so inho- 
mogeneities will not have had time to phase mix and con- 
tinued disturbance from mergers will not have relaxed (see 
Tremaine 1992 for additional discussion). It is likely that in- 
trinsic fluctuations may play a significant role in long-term 
galaxian evolution. In particular, a companion paper (Wein- 
berg 1997) elaborates the nature of satellite halo excitation 
and its role in producing warps using a similar formalism. N- 
body simulation of this process revealed that intrinsic noise 
had a similar effect on the disk. It is difficult to deduce 
the power produced by phase mixing streams based on the 
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point-mass noise from 10 s , 10 6 or 10 7 Mq fuzzy blobs and 
the degree of granularity will depend on the nature of the 
dark matter and halo formation history. Speculating further 
nonetheless, it is conceivable that a wide variety of observed 
disk structure such as arms and bars may have its roots in 
halo structure. 
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